4, 2008 16:53 Proceedings Trim Size: 9in x 6in 



cardall'Openlssues 



THE LONG TERM: SIX-DIMENSIONAL CORE-COLLAPSE 
SUPERNOVA MODELS 



C. Y. CARDALL 1 ' 2 , A. O. RAZOUMOV 1 ' 2 ' 3 , E. ENDEVE 1 ' 2 ' 3 , AND 
A. MEZZACAPPA 1 

1 Physics Division, 
Oak Ridge National Laboratory, 
Oak Ridge, TN 37831-6354, USA 

2 Department of Physics and Astronomy, University of Tennessee, Knoxville, 

TN 37996-1200, USA 

3 Joint Institute for Heavy Ion Research, 
Oak Ridge National Laboratory, 
Oak Ridge, TN 37831-6374, USA 

The computational difficulty of six-dimensional neutrino radiation hydrodynamics 
has spawned a variety of approximations, provoking a long history of uncertainty in 
the core-collapse supernova explosion mechanism. Under the auspices of the Teras- 
cale Supernova Initiative, we are honoring the physical complexity of supernovae 
by meeting the computational challenge head-on, undertaking the development 
of a new adaptive mesh refinement code for self-gravitating, six-dimensional neu- 
trino radiation magnetohydrodynamics. This code — called GenASiS, for General 
Astrophysical Simulation System — is designed for modularity and extensibility of 
the physics. Presently in use or under development are capabilities for Newto- 
nian self-gravity, Newtonian and special relativistic magnetohydrodynamics (with 
'realistic' equation of state), and special relativistic energy- and angle-dependent 
neutrino transport — including full treatment of the energy and angle dependence 
of scattering and pair interactions. 



1. The Challenges of Core-collapse Supernovae 

In taking stock of 'long-term' efforts to understand core-collapse super- 
novae, we reflect upon the fact that supernovae have been challenging us 
for centuries. Their very existence helped overturn worldviews. Their ex- 
plosion mechanisms and remnants involve all four fundamental forces, and 
many (if not most) branches of physics. A cornucopia of electromagnetic ra- 
diation observables continues to provide intriguing puzzles, and yields some 
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clues regarding the violent proceedings of a massive star's death. But direct 
observational penetration of the secrets of neutron star birth and initiation 
of the explosive ejection of the stellar envelope demands extraordinary ef- 
forts aimed at the detection of gravitational waves and neutrinos, the only 
messengers carrying direct information from the extreme conditions of a 
newly-collapsed stellar core. And the associated theoretical penetration — 
required for both the prediction and interpretation of expected gravitational 
wave and neutrino signals — comprises algorithmic and computational issues 
that will challenge computational physicists and tax state-of-the-art super- 
computers for years to come. 

In western civilization, supernovae played a role in changing prevailing 
notions of the universe in at least two eras. Remarkably, of the handful 
of supernovae in our Milky Way Galaxy recorded by humanity, two were 
observed by Tycho (1572) and Kepler (1604). Tycho's detailed observations 
established that the 'new and never previously seen star' of 1572 — and also 
another transient celestial phenomenon, a comet of 1577 — were beyond the 
moon's orbit, 'new phenomena in the ethereal world,' contributing to the 
overthrow of the Aristotelian worldview that included immutable heavens. 
In modern times, supernovae figured in the debate over whether the spiral 
nebulae were separate galaxies, each an 'island universe' comparable to 
our Milky Way. It was recognized that the 'novae' or 'new stars' seen in 
these nebulae would have to be much more luminous than typical novae 
occuring in our galaxy. The phrases "giant novae," novae of "impossibly 
great absolute magnitudes," "exceptional novae," and the German term 
"Hauptnovae" or "chief novae" were used during the 1920s. 1 In a review 
article, Zwicky explained that it was deduced that 'supernovae' were about 
a thousand times as luminous as 'common novae,' and claimed that "Baade 
and I first introduced the term 'supernovae' in seminars and in a lecture 
course on astrophysics at the California Institute of Technology in 1931." 2 

Supernovae are classified by astronomers into two broad classes based 
on their optical spectra. 3 These classes are 'Type I,' which have no hydro- 
gen features, and 'Type II,' which have obvious hydrogen features. These 
types have further subcategories, depending on the presence or absence of 
silicon and helium features in Type I, and the presence or absence of nar- 
row hydrogen features in the case of Type II. In particular, supernovae of 
Type la exhibit strong silicon lines, those of Type lb have helium lines, and 
those of Type Ic do not have either of these. Astronomers have also iden- 
tified a number of distinct characteristics in supernova light curves (total 
luminosity as a function of time) . 
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There are two basic physical mechanisms for supernovae, but these do 
not line up cleanly with the observational categories of Type I and Type II. 
Type la supernovae are caused by a thermonuclear runaway that consumes 
an entire white dwarf, thought to be induced by accretion of matter from a 
companion star. Supernovae of Type lb, Ic, and II are produced by a totally 
different mechanism: the catastropic collapse of the core of a massive star. 
The observational distinctions of presence or absence of hydrogen or helium 
turn out to be unrelated to the mechanism; they depend on whether the 
outer hydrogen and helium layers of the star — which have nothing to do 
with the collapsing core — have been lost to winds or accretion onto a binary 
companion during stellar evolution. Of the two physical mechanisms, core- 
collapse supernovae are the focus of the present discussion. 

We now consider the core-collapse supernova process in more detail. 
Shortly after the discovery of the neutron in the early 1930s, Baade and 
Zwicky declared, "With all reserve we advance the view that supernovae 
represent the transitions from ordinary stars to neutron stars, which in their 
final stages consist of extremely closely packed neutrons." 4 This turned 
out to be true, at least for some 'core-collapse' supernovae (those of Type 
Ib/Ic/II); a black hole is another possible outcome. The dominant flcshing- 
out of the core collapse process in the last two decades 3, has been the delayed 
neutrino-driven explosion mechanism. 6,7 

A core-collapse supernova results from the evolution of a massive star. 
For most of their existence, stars burn hydrogen into helium. In stars at 
least eight times as massive as the Sun (8 M ), temperatures and densi- 
ties become sufficiently high to burn through carbon to oxygen, neon, and 
magnesium; in stars of at least ~ 10 Mq , burning continues through silicon 
to iron group elements. The iron group nuclei are the most tightly bound, 
and here burning in the core ceases. 

The iron core — supported by electron degeneracy pressure instead of 
gas thermal pressure, because of cooling by neutrino emission from carbon 
burning onwards — eventually becomes unstable. Its inner portion under- 
goes homologous collapse (velocity proportional to radius), and the outer 
portion collapses supersonically. Electron capture on nuclei is one insta- 
bility leading to collapse, and this process continues throughout collapse, 
producing neutrinos. These neutrinos escape freely until densities in the 
collapsing core become so high that even neutrinos are trapped. 

Collapse is halted soon after the matter exceeds nuclear density; at 



a See for example Ref. 5 for some information on earlier views of the mechanism. 
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this point ("bounce"), a shock wave forms at the boundary between the 
homologous and supersonically collapsing regions. The shock begins to 
move out, but after the shock passes some distance beyond the surface of 
the newly-born neutron star, it stalls as energy is lost to neutrino emission 
and endothcrmic dissociation of heavy nuclei falling through the shock. 

It is natural to consider neutrino heating as a mechanism for shock 
revival, because neutrinos dominate the energetics of the post-bounce evo- 
lution. Initially, the nascent neutron star is a hot thermal bath of dense 
nuclear matter, electron/positron pairs, photons, and neutrinos, contain- 
ing most of the gravitational potential energy released during core collapse. 
Neutrinos, having the weakest interactions, are the most efficient means of 
cooling; they diffuse outward on a time scale of seconds, and eventually 
escape with about 99% of the released gravitational energy. 

Because neutrinos dominate the energetics of the system, a detailed un- 
derstanding of their evolution will be integral to definitive accounts of the 
supernova process. If we want to understand the origin of the explosion 
with energy ~ 10 51 erg, we cannot afford to lose (or gain) more than this 
amount during the period covered by the simulation. This requires careful 
accounting of the neutrinos' much larger contribution to the system's en- 
ergy budget. (For further discussion, and a review of work recognizing the 
importance of this point, see Ref. 8 ). 

What sort of computation is needed to follow the neutrinos' evolution? 
Deep inside the newly-born neutron star, the neutrinos and the fluid are 
tightly coupled (nearly in equilibrium); but as neutrinos are transported 
from inside the neutron star, they go from a nearly isotropic diffusive regime 
to strongly forward-peaked free-streaming. Heating behind the shock oc- 
curs precisely in this transition region, and modeling this process accurately 
requires tracking both the energy and angle dependence of the neutrino dis- 
tribution functions at every point in space. 

A full treatment of this six-dimensional neutrino radiation hydrodynam- 
ics problem is a major challenge, too costly for contemporary computational 
resources. While much has been learned over the years through simulation 
of model systems of reduced dimensionality, there is as yet no robust con- 
firmation of the delayed neutrino-driven scenario described above (see Sec. 

Recent detections of a handful of unusually energetic Type Ib/c su- 
pcrnovae (often called 'hypernovae') in connection with gamma-ray bursts 
pose additional challenges to theory and observation. Prominent exam- 
ples of this supernova/gamma-ray burst connection include SN1998bw / 
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GRB980425, 9 * 10 SN20021t / GRB021211, 11 SN2003dh / GRB030329, 12 * 13 
and SN20031w / GRB031203; 14 ' 15 ' 16 ' 17 there arc probably many others (see, 
for example, Ref. 18 ). Like many gamma-ray bursts without direct evidence 
for a supernova connection, GRB030329 has evidence of a jet; GRB980425 
and GRB031203 do not, and are also undcrluminous gamma-ray bursts 
(but still unusually energetic Type Ib/c supernovae). 19 ' 20 Determining the 
relative rates of jet-like hypernovae, non-jet-like hypernovae, and 'normal' 
supernovae — and the possible associated variety of mechanisms — are im- 
portant challenges. 

In summary, the details of how the stalled shock is revived sufficiently to 
continue plowing through the outer layers of the progenitor star are unclear. 
In normal supernovae, it may well be that some combination of neutrino 
heating of material behind the shock, convection, and instability of the 
spherical accretion shock leads to the explosion (see Sec. 2). It is tempting 
to think that rotation (for example, Refs. 21 > 22 ) and magnetic fields (for 
example, Ref. 23 ) in more massive progenitors may play a more significant 
role in the rare jet-like hypernovae, perhaps giving birth to 'magnetars,' the 
class of neutron stars with unusually large magnetic fields. This temptation 
appears to be sweetened by observational support. 24 ' 25 (Observations of two 
nearby supernova remnants may suggest that rotation and magnetic fields 
also operate in normal supernovae, perhaps subdominantly. 26 ) 

From the above discussion, several key aspects of physics that a core- 
collapse simulation must address can be identified; these are discussed in 
sections that follow, after a discussion of the history of approximate treat- 
ments of neutrino radiation transport and an overview of our new code. 

2. History of Neutrino Radiation Hydrodynamics 

While in general terms supernovae have been challenging us for centuries, 
the challenge of their simulation via computer modeling has 'only' been 
with us for a few decades — a 'short term' in comparison with centuries, but 
still a 'long term' in comparison with the time scales of individual academic 
careers. 

Here we sketch the last two decades' progress on one critical aspect of 
core-collapse supernova simulations: the high dimensionality (three space 
and three momentum space dimensions — not to mention time dependence) 
of neutrino radiation hydrodynamics (sec Table 1). The development of 
this aspect of the simulations is intertwined with important advances in 
the field, but of course does not represent every insight relevant to the 
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explosion mechanism obtained via simulation or otherwise. 

We pick up the story in 1982, when simulations showing the stalled shock 
reenergized by neutrino heating on a time scale of hundreds of milliseconds 
were first performed. 6 This was initially achieved in a simulation with a 
total of 2 dimensions (spherical symmetry and energy-dependent neutrino 
transport). But these simulations required significant rezoning, possibly at- 
tended by nontrivial numerical error; 6 and further, with the introduction of 
full general relativity and a correction in an outer boundary condition, 34 it 
became clear that these models would not explode without a mock-up of a 
doubly-diffusive fluid instability in the newly-born neutron star that serves 
to boost neutrino luminosities 34 ' 35 ' 36 ' 37 — a simulation of effective total di- 
mensionality "2.5" (see Table 1). That the necessary conditions exist for 
this particular instability to operate has been disputed; 41 ' 62 ' 63 and though 
related phenomena may operate, 63 more recent simulations with energy- 
dependent neutrino transport and true two-dimensional fluid dynamics in- 
dicate that fluid motions are either suppressed by neutrino transport 42 or 
have little effect on neutrino luminosities and supernova dynamics. 59 ' 61 

Recognizing that the profiles obtained in spherically symmetric simu- 
lations implied convective instabilities, and that observations of supernova 
1987A also pointed to asphericities, several groups explored fluid motions 
in two spatial dimensions in the supernova environment in the 1990s. In 
two spatial dimensions, the computational limitations of that era required 
approximations that simplified the neutrino transport. 

One class of simplifications allowed for neutrino transport in "1.5" or 
2 spatial dimensions, but with neutrino energy and angle dependence inte- 
grated out, reducing a five dimensional problem to "1.75" or 2 effective to- 
tal dimensions (see Table l). 37 > 38 . 39 These simulations exhibited explosions, 
and elucidated an undeniably important physical effect: a negative entropy 
gradient behind the stalled shock results in convection that increases the 
efficiency of heating by neutrinos. However, in the scheme of Table 1, the 
inability to track the neutrino energy dependence in these simulations could 
be viewed as a minor step backwards in effective total dimensionality. The 
energy dependence of neutrino interactions has the important effect of en- 
hancing core deleptonization, which makes explosions more difficult; 30 ' 64,65 
this raised the question of whether the exploding models of the early- and 
mid-1990s were too optimistic. 

This concern about the lack of neutrino energy dependence received 
some support from a simulation in the late 1990s involving a different simpli- 
fication of neutrino transport: the imposition of energy-dependent neutrino 



Tabic 1. Selected neutrino radiation hydrodynamics milestones in stellar collapse simulations studying the long-term fate of the shock. 



Group 


Year 


Explosion 


lotai 


Fluid space 


v space 


v momentum 








dimensions 


dimensions 


dimensions 


space dimensions 


Lawrence Livcrmorc 2 '' 6 '' 


1982 


Yes* 


2 


1 (PN) 


1 


1 {0(v/c)) 


Lawrence Livcrmorc 28,29 


1985 


Yes* 


"2.25" 


"1.5" NS (PN) 


1 


1 (0(v/c)) 


Florida Atlantic 30 - 31 - 32 - 33 


1987 


No 


2 


1 (GR) 


1 


1 (0(v/c)) 


Lawrence Livcrmorc 34 - 35 ' 30 


1989 


Yes* 


"2.25" 


"1.5" NS+HR (GR) 


1 


1 (GR) 


Lawrence Livermore 3 ' 


1992 


Yes* 


2 


2 HR (N) 


2 


(N) 


Los Alamos 38 


1993 


Yes* 


"1.75" 


2 (N) 


"1.5" thick/thin 


(PN) 


Arizona 39 


1994 


Yes* 


"1.75" 


2 (N) 


"1.5" ray-by-ray 


(N) 


Florida Atlantic 4 "- 41 


1994 


No 


"2.25" 


"1.5" NS (GR) 




1 (0(v/c)) 


Oak Ridge 42 - 43 


1996 


No 


"2.5" 


2 (N) 




1 (0(v/c)) 


Max Planck 44 ' 45 ' 4 "' 47 


2000 


No, Yes* (ONeMg) 


3 


1 (N) 




2 (0(v/c)) 


Oak Ridge 48 ' 49 ' 50 ' 51 ' 52 


2000 


No 


3 


1 (N) 




2 (0(v/c)) 


Arizona^ ! 


2002 


No 


3 


1 (N) 




2 (O(v/o)) 


Oak Ri dge 4«,4«,50,55,5ti,57 


2000 


No 


3 


1 (GR) 




2 (GR) 


Los Alamos 38 - 58 


2002 


Yes* 


"2.5" 


3 (N) 


"2" thick/thin 


(PN) 


Max Planck 45 ' 5 ^' til ' 4V 


2003 


No, Yes* (180°) 


"3.75" 


2 (PN) 


"1.5" ray-by-ray 


2 (0(v/c), PN) 



Note: The "Yes" entries in the "Explosion" column are all marked with an asterisk as a reminder that questions about the simulations — 
described in the main text — have prevented a consensus about the explosion mechanism. "Total dimensions" is the average of "Fluid 
space dimensions" and "v space dimensions," added to "v momentum dimensions." The abbreviation "N" stands for 'Newtonian,' while 
"PN" — for 'Post-Newtonian' — stands for some attempt at inclusion of general relativistic effects, and "GR" denotes full relativity. A 
space dimensionality in quotes — like "1.5" — denotes an attempt at modeling higher dimensional effects within the context of a lower 
dimensional simulation. For the fluid, this is a mixing- length prescription in the neutron star ("NS") or the heating region ("HR") 
behind the stalled shock. For neutrino transport, it indicates one of two approaches: multidimensional diffusion in regions with strong 
radiation/fluid coupling, matched with a spherically symmetric 'light bulb' approximation in weakly coupled regions ("thick/thin"); or 
the (mostly) independent application of a spherically symmetric formalism/algorithm to separate spatial angle bins ( "ray- by-ray" ) . 
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distributions from spherically symmetric simulations onto fluid dynamics 
in two spatial dimensions. 43 Unlike the simulations discussed above, these 
did not explode, casting doubt upon claims that convection-aided neutrino 
heating constituted a robust explosion mechanism. 

The nagging qualitative difference between spatially multidimensional 
simulations with different neutrino transport approximations motivated in- 
terest in the possible importance of even more complete neutrino trans- 
port: Might the retention of both the energy and angle dependence of 
the neutrino distributions improve the chances of explosion, as preliminary 
"snapshot" studies suggested? 66 ' 53 Of necessity, the first such simulations 
were performed in spherical symmetry, which nevertheless represented an 
advance to a total dimensionality of 3 (see Table 1). Results from three 
different groups are in accord: Spherically symmetric models of iron core 
collapse do not explode, even with solid neutrino transport 44 ' 52,54 and gen- 
eral relativity. 55 ' 57 Recently, however, it has been shown that the more 
modest oxygen/neon/magnesium cores of the lightest stars to undergo core 
collapse (8-10 M ) may explode in spherical symmetry. 46 ' 47 

The current state of the art in neutrino transport in supernova simu- 
lations determining the long-term fate of the shock has been achieved by 
a group centered at the Max Planck Institute for Astrophysics in Garch- 
ing, who deployed their spherically symmetric energy- and angle-dependent 
neutrino transport capability 45 along separate radial rays, with partial cou- 
pling between rays. 60 Initial results — from axisymmetric simulations with 
a restricted angular domain — were negative with regards to explosions (in 
spite of the salutary effects of convection, and also rotation), 59 apparently 
supporting the results of Ref. 43 . An explosion was seen in one simulation 67 
in which certain terms in the neutrino transport equation corresponding to 
Dopplcr shifts and angular aberration due to fluid motion were dropped; 
this simulation also yielded a neutron star mass and nucleosynthetic con- 
sequences in better agreement with observations than the "successful" ex- 
plosion simulations of the 1990s, 38 ' 39 arguably because of more accurate 
neutrino transport in the case of both observables. The continuing lesson 
is that getting the details of the neutrino transport right makes a difference. 

In addition to accurate neutrino transport, low-mode (£ = 1,2) insta- 
bilities that can develop only in simulations allowing the full range of polar 
angles may make a subtle but decisive difference, as in an explosion recently 
reported by the Garching group. 61,47 This achievement was presaged by 
earlier studies of the supernova context, which featured a demonstration 
of the tendency for convective cells to merge to the lowest order allowed 
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by the spatial domain 68 and a ncwly-rccognizcd spherical accretion shock 
instability 69 (discovered independently in a different context in Ref. 70 ). 
These these global asymmetries may be sufficient to account for observed 
asphericities that have often been attributed to rotation and/or magnetic 
fields. 

Surely every 'Yes' entry in the explosion column of Table 1 has been 
hailed in its time as 'the answer' (at least by some!), and as a community 
we cannot help hoping once again that these recent developments mark the 
turning of a corner; but important work remains to verify if this is the case. 
Several groups are committed to further efforts. For example, the Terascale 
Supernova Initiative (TSI, which includes authors of Ref. 69 ) comprises ef- 
forts aimed at 'ray-by-ray' simulations 71 like those of the Garching group, 
as well as full spatially multdimensional neutrino transport, both with en- 
ergy dependence only 72 and with energy and angle dependence (Sec. 6, and 
Ref 73 ). Delineation of the possible roles of rotation and magnetic fields 
are also being pursued by TSI. At least one other group is pursuing full 
spatially multidimensional neutrino transport. 74 ' 75 

3. GenASiS: Philosophy and Basic Features 

As discussed in Sec. 1, a core-collapse supernova involves a six-dimensional 
radiation hydrodynamics problem, making it a major computational chal- 
lenge. Even three-dimensional pure hydrodynamics problems (with only 
space dimensions, no momentum space) have only become relatively com- 
mon in the last few years, with manageable workflows on today's terascale 
machines (<~ 10 12 bytes of memory and flop/s). To begin to get a feel for the 
requirements of radiation hydrodynamics, consider just a five-dimensional 
problem, in which axisymmctry in the space dimensions is assumed. For 
example, supposing the numbers of spatial zones in spherical coordinates 
(r,0) to be (256, 128), and the numbers of momentum bins in energy and 
angle variables (e,i9,ip) to be (64,32,16), of order 10 10 bytes are required 
just to store one copy of one neutrino distribution function. While this 
gives rise to a taxing (but not necessarily insurmountable) workflow on 
terascale machines, it is apparent that the addition of the third spatial 
dimension — necessary for full exploration of the interacting effects of con- 
vection, rotation, and magnetic fields — will require petascale systems. But 
petascale systems will eventually be available (five to seven years is the 
current expectation); and given the long development time scales of so- 
phisticated software, we believe it wise to develop our code with the full 
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six-dimensional capability, even if it is only deployed in five dimensions in 
the near term. 

The computational demands of radiation hydrodynamics can be ame- 
liorated by "adapative mesh refinement" (AMR). The basic idea of AMR is 
to employ high resolution only where needed in order to conserve memory 
and computational effort. Our current expectation is to allow for refine- 
ment only in the space dimensions. This will help with the management of 
two difficulties: the large dynamic range in length scales associated with the 
density increase of six orders of magnitude that occurs during core collapse, 
and adequate resolution of particular features of the flow (the shock, for 
instance). With Eulerian codes in multiple space dimensions, these tasks 
require high resolution; and particularly for radiation hydrodynamics, the 
savings achievable by reducing the number of zones is considerable, since an 
entire three-dimensional momentum space is carried by each spatial zone. 

Of the two basic types of AMR on structured grids — the block- 
structured and zone-by-zone varieties — we have chosen the zone-by-zone 
approach for use with neutrino radiation hydrodynamics. Block-structured 
AMR 76 involves the deployment of subgrids of a certain reasonable mini- 
mum size (e.g. eight or sixteen zones per side) at various levels of refine- 
ment. A basic solver routine is applied independently to each subgrid. Ex- 
tra spatial zones (referred to as 'guard zones') are required on the edges of 
each subgrid, which carry information from neighboring subgrids; these be- 
come the boundary conditions applied by the solver routine. The strategy is 
designed for explicit solution algorithms, in which the functions describing 
time evolution need only be evaluated at the previous time step. However, 
the rapid time scales of neutrino interactions with the fluid require an im- 
plicit solution algorithm, in which the functions describing the evolution 
of the neutrino radiation field are evaluated at the current (that is, new) 
time step. Because of this mismatch with the intended purposes of block- 
structured AMR (implicit vs. explicit evolution), and the fact that popular 
block-structured AMR community packages did not seem readily amenable 
to handling momentum space variables in a natural way, we decided upon 
another flavor of AMR: the zone- by-zone refinement approach. 77 In this 
method, individual zones are refined (typically by bisection) and coarsened 
as needed. This provides more flexibility than the block-structured ap- 
proach; the fine-grained control allows for maximum savings in the number 
of spatial zones deployed. A drawback for many users is that a single-grid 
explicit solver cannot be used "as is." Instead, new solution algorithms 
must be developed that address the entire hierarchical data structure (Rcf. 
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77 and our Sec. 4 are examples for explicit hydrodynamics); but we are 
required to develop such 'global solvers' for gravity (Sec. 5) and implicit 
neutrino transport (Sec. 6) anyway. And as a bonus, the need to carry 
memory-wasting 'guard zones' is obviated. 

In implementing the zone-by-zonc-refinement approach to the represen- 
tation of spacetime, we have tried to follow object-oriented design principles 
to the extent allowed by Fortran 90/95. Figure 1 outlines the basic data 
structures we use to model the ideal of a continuous spacelike slice with a 
discretized approximation. (The hierarchy of structures, and our operations 
on them with well-controlled interfaces, are instances of the object-oriented 
principles of inheritance and encapsulation.) A region of a spacelike slice 
is represented by an object of zone Array Type. Each such object contains 
an array of objects of zoneType, along with information about the coordi- 
nates of the zones and pointers to neighboring zone arrays. Each zone, an 
object of zoneType, contains various forms of stress-energy, each of which 
is a separate object. Figure 1 shows a perfect fluid and a radiation field; in 
the code we have an electromagnetic field as well. Each zone has a pointer 
to another object of zoneArrayType, whose allocation constitutes refine- 
ment of that zone; this structure can be extended to arbitrarily deep. A 
simple two-dimensional pure hydrodynamics test problem computed with 
our adaptive mesh code is shown in Fig. 2. 

The word "General" that goes into the name of our code, GenASiS — for 
General ^Istrophysical Simulation System — may give an initial impression 
of a messianic quest to create an impossibly all-purpose code for solving all 
conceivable problems in astrophysics and cosmology; but the code's 'gen- 
erality' is, of course, considerably more modest: It refers to the use of 
Fortran 90's facility for function overloading. (This is an instance of the 
object-oriented principle of polymorphism.) This allows a generic function 
name to have several different implementations, providing for extensibility 
of the physics: Different equations of state, hydrodynamic flux methods, co- 
ordinate systems, gravity theories, and so forth can be employed by adding 
new implementations of generic function names, without having to go back 
and change basic parts of the code to implement new physics. 

4. Magnetohydrodynamics 

We employ a conservative formulation of the equations of magnetohydro- 
dynamics. The Newtonian case will be described here. Conservation of 
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1 Radiation 


1 ZoneArravType | 
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ZoneTvpe 


PerfectFluid 


II PerfectFluid II 










1 Radiation 


1 Radiation 


1 ZoneArravType 1 


1 ZoneArravType 







Figure 1. Data structures used in an adaptive mesh for radiation hydrodynamics. 




Figure 2. Density in a two-dimensional generalization of the shock tube. Red and blue 
indicate high and low density respectively. Left: Initial state. Right: Evolved state. 
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baryons is described by the equation 

dn _ , 

— + V • (nv) = 0, 

where n is the baryon number density and v is the fluid velocity, 
equation 



9 < ^ „ 
— (mnv) + V 



ran vv + I p 



BB 



-mnV$ 



(1) 
The 

(2) 



describes conservation of momentum. Here m is the average baryon mass, 
p is the pressure, B is the magnetic field, $ is the gravitational potential 
(discussed further in Sec. 5), and 1 is the unit tensor. Units of the magnetic 
field arc chosen such that the vacuum magnetic permeability is unity. One 
way to express conservation of energy is 



d_ 

di 



mn 



2 



(e + p + B 2 )v + ™ (v 2 + $) v - B(v ■ B) 



mn 

dm 
~dt 



V * + v V$) - 
v ■ V/// | . (3) 



where e is the internal energy density, and * is a kind of "gravitational 
vector potential" (also discussed in Sec. 5). We have opted to employ the 
baryon number density, instead of the mass density, in explicit deference 
to the fact that mass is not conserved in the presence of nuclear reactions. 
The energy input from nuclear reactions has then resulted naturally in the 
derivation of Eq. (3), appearing in the last two terms of the right-hand 
side. We hope to add nuclear reaction networks to our code in the future. 

This formulation is called "conservative" because volume integrals of 
the divergences in Eqs. (l)-(3) are related to surface integrals through the 
divergence theorem: 



dV(V -F)= d> F dA. (4) 

'V JdV 

The physical meaning of a conservative equation is that (modulo source 
terms) the time rate of change of a conserved quantity in a volume is equal 
to a flux F through the volume's enclosing surface. This meaning is built 
into the finite-difference representation of Eqs. (l)-(3); divergences are 
represented in discrete correspondence to their mathematical definition, 
using zone volumes V and face areas A: 



v • F -+ vb E [( A » F v - ^ pq ^ 



(5) 
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Here q runs over the three space dimensions. A double-headed arrow (<-») 
indicates evaluation at a zone center. Left-arrows (<— q) and right-arrows 
(q — >) denote evaluation at zone inner and outer faces respectively, in the 
q direction; dimensions other than q are evaluated at the zone centers. 
In this way the divergence theorem is replicated in every zone, ensuring 
global conservation to machine precision. Use of generalized zone volumes 
and areas in Eq. (5) enables the use of curvilinear coordinates ("ficticious 
forces" arising from curvilinear coordinates must also be included in the 
momentum equation). 

The evolution of the magnetic field is described by Faraday's law: 



supplemented by the constraint 



In Eq. (6), we take the electric field to be E = —v x B, in accordance with 
the usual astrophysical assumption of a perfectly conducting medium. 

While Eq. (6) for the evolution of the magnetic field does not require 
conservation of the magnetic field, Eq. (7) requires the magnetic field 
to be divergence-free at all times. In the presence of discontinuous flow 
numerical solutions to Eqs. (l)-(6) can produce severe unphysical artifacts 
if this requirement is not met, 78 but it can be automatically enforced by 
the method of constrained transport. 79 Integrating Eq. (6) over a zone's 
enclosing surface, the left-hand side becomes the volume integral of V • B, 
via the divergence theorem. The right-hand side is a sum over area integrals 
over each zone face. With Stokes' theorem, 



each of these surface integrals becomes a line integral around the zone face 
boundary. Summed over all faces, two line integrals in opposite directions 
cancel on every zone edge, enforcing V • B = in the zone as desired. The 
method of constrained transport, then, is to evaluate V x E on zone edges 
in discrete correspondence to the mathematical definition of the curl, using 
zone face areas A and edge lengths L: 



V 



B = 0. 



(7) 




(8) 



(9) 
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Here r runs over the two space dimensions orthogonal to a particular direc- 
tion q, and s indicates the direction perpendicular to both q and r. Left- 
arrows (<— r) and right-arrows (r — ►) denote evaluation along face inner 
and outer edges respectively, in the r direction. This ensures divergence- 
free evolution of the discrete representation of the area-averaged magnetic 
field, with components located on the appropriate zone faces, for all times 
to machine precision, provided the initial magnetic field satisfies Eq. (7). 

Accurate computation of fluxes at zone faces and electric fields at zone 
edges is a key feature. So-called "central schemes" have been noted recently 
by astrophysicists for their ability to capture shocks with an accuracy com- 
parable to Riemann solvers, but with much greater simplicity. 80 In partic- 
ular, we employ so-called "HLL" versions of these schemes for both fluid 
conservation laws 81 and the magnetic induction equation. 82 We achieve 
second order in space by linear interpolation within zones (as usual, a 
slope limiter — deployed where necessary in order maintain discontinuities — 
reduces the treatment to first order). 

While we have taken Khokhlov's zone-by-zone refinement approach 77 as 
our basic paradigm, we have made a novel extension to evolve the magnetic 
field, and use a different time-stepping scheme. As in Khokhlov's work, 
fluxes are computed at each zone interface only once, with the results used 
to update zones on both sides of the interface. At coarse/fine interfaces, 
fluxes are computed only on the faces of the refined zones. We have de- 
veloped a similar approach for the induction equation: the electromotive 
forces on the zone edges are computed only once, and are used to update all 
zones sharing that edge. We use a second-order Runge-Kutta time stepping 
algorithm, made possible by the semi-discrete formulation of the central 
scheme. In doing so we evolve all levels of the mesh synchronously, unlike 
Khokhlov's approach of evolving refined levels with greater frequency. In 
addition to making it possible to take advantage of the semi-discrete for- 
mulation for time evolution, problems we encountered in self-gravitating 
systems with 'asynchronous' evolution a la Khokhlov were avoided. 

Parallelization is achieved by giving each processor its share of spatial 
zones. Partitioning is accomplished by walking through all levels of the 
mesh in a recursive manner similar to a Morton space-filling curve; the 
result is a mapping of the multidimensional mesh to a one dimensional 
"string" of zones, which, when cut into pieces of uniform length, leaves 
each processor with roughly the same number of zones at each level of 
refinement. 

An equation of state determines p, the quantity in Eqs. (l)-(3) whose 
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determination has not yet been mentioned. So far, the code has two over- 
loaded options. One is the familiar polytropic equation of state: 

p = o r , (10) 

e = (r-i)"V (ii) 

where T is a specified parameter, and k is updated in response to changes 
in e determined from Eq. (3). We have also implemented a "realistic" 
equation of state 83 suitable for problems involving nuclear matter. This 
equation of state takes as input the temperature, baryon number density, 
and electron fraction Y e , defined by 

Y e = I!£__^± (12) 
n 

where n e - and n e + are the number densities of electrons and positrons 
respectively. When using this "realistic" equation of state, an advection 
equation for Y e must be added to the above list of conservation laws. 

We have been working with a number of hydrodynanic and magne- 
tohydrodynamic test problems, one of which is shown here. The rotor 
problem — which consists of a rapidly rotating dense fluid, initially cylindri- 
cal, threaded by an initially uniform magnetic field — was devised to test the 
onset and propagation of strong torsional Alfven waves into the ambient 
fluid. We have computed a version of the rotor problem with initial data 
identical to a so-called 'second rotor problem,' 85 and display the results in 
Fig. 3. 



5. Newtonian Gravity 

The Poisson equation for the Newtonian gravitational potential $ is 

V 2 $ = 47rGmn, (13) 

where G is the gravitational constant. Our finite-differenced approach to 
Eq. (13) is based on the fact that the Laplacian is a divergence of a gradient: 

V • V$ = 4ttG ran. (14) 

We discrctizc this equation in a manner similar to Eq. (5). There results 
a linear system for the values of $ at the center of every leaf zone. In the 
matrix representation of this linear system, each row corresponds to the 
discrete version of Eq. (14) centered on a given "leaf zone" (those not having 
refined children). On a single- level grid, the resulting matrix would have 
three, five, and seven bands in one, two, and three dimensions respectively. 
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Figure 3. Density (left panel) and thermal pressure (right panel) at t = 0.295 for the 
second rotor problem given in Ref. 85 . 40 countours were used to produce the plots, 
with 0.512 < ran < 9.622 and 0.010 < p < 0.776. A 200 X 200 grid was used to produce 
the results. 



With adaptive mesh refinement, the matrix structure becomes more diffuse, 
because several refined zones may contribute to the discrete representation 
of V<f> at a zone face featuring a coarse/fine interface. Each processor fills 
in the portion of the matrix corresponding to its share of zones, and we rely 
on the PETSc library (http://www-unix.mcs.aiil.gov/petsc/petsc-2/) 
to perform the distributed sparse matrix inversion. 

We now discuss the contribution of gravity to the fluid energy evolution 
in Eq. (3). Without the gravitational potential energy included in the 
"total energy" in the time derivative, Eq. (3) appears in the more familiar 
form 



d_ 

dt 



e + -mnv x 



+ V- 



(e + p + -mnv 2 )v 



-mnv ■ V$ 



Including the graviational potential energy density (mn/2)$ in the time 
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derivative, and making use of baryon conservation, we have 
d_ 

at 

,,2 



V • 



(e +p)v 



ran 
ran 



(v 2 + $) 



+ 



v- 



mn ( 9$ 



- v ■ V$ - 



dt 

( dm 

— v ■ Vto 

V dt 



Using the formal solution for 

*(x,t) = -G J 



mix 1 ',t)n(x' \t) d 3 x' 



(16) 



(17) 



\x — x' 

together with baryon conservation (and neglecting time derivatives of to), 
one can show that 

~dt 

where satisfies a vector Poisson equation: 

V 2 * = AirGmnv. 



(18) 



(19) 

This may be solved in a manner similar to that used to solve for $ (though 
the vector Poisson equation contains some additional terms in curvilinear 
coordinates). Conservation of total energy, including gravitational, is not 
local: the gravitational source terms (the first two terms on the right-hand 
side of Eq. (3)) vanish only upon integration over all space. 



6. Neutrino Radiation Transport 

Here we briefly describe our approach to the greatest computational chal- 
lenge in supernova simulations: neutrino radiation transport. Neutrino 
distributions must be tracked in order to compute the transfer of lepton 
number and energy between the neutrinos and the fluid. There are three 
major challenges. One challenge is constructing a discretization that allows 
both energy and lepton number to be conserved to high precision. The two 
other challenges are associated with the limits of computational resources: 
the solution of a very large nonlinear system of equations, and neutrino 
interaction kernels of high dimensionality 

Energy conservation is an obvious measure of quality control, and care 
with the transport formalism and differencing can help achieve it. The 
importance of energy conservation is brought into focus by this question: 
How should we interpret the prediction of a ~ 10 51 erg explosion in a model 
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where the total energy varies during the course of the simulation by ~ 10 51 
erg or more? To achieve the required precision (say global energy changes of 
less than ~ 10 50 erg), conservative formulations are a useful starting point, 
and relativistic treatments avoid quantitatively non-negligible conflicts at 
0(v 2 /c 2 ) between the number and energy transport equations. 86 > 57 > 87 > 8 
Finite-differencing that simultaneously satisfies energy and lepton number 
conservation has been implemented in spherical symmetry, 57 and should be 
pursued in multiple spatial dimensions as well. 

Our algorithm for solving the large nonlinear system has been described 
elsewhere. 88 ' 73 The large system of equations requiring inversion (as op- 
posed to explicit updates) results from the disparity between hydrodynamic 
and particle interaction time scales, which motivates implicit time evolu- 
tion. The nonlinear solve is achieved with the Newton-Raphson method. 
A fixed-point method employing a preconditioner that splits the space and 
momentum space couplings is used for the linear solve required within each 
Newton-Raphson iteration. 88 An advantage of this linear solver method is 
that the dense blocks representing couplings in momentum space — which 
cannot all be stored at once — need only be constructed a few at a time, 
used in all steps required in a given fixed-point iteration, and discarded. In 
contrast, other linear solver algorithms seem to require dense blocks to be 
discarded and rebuilt multiple times in each iteration. We have sucessfully 
tested this solver on a two-dimensional problem in spherical coordinates 
with a static background and a simple emission/ absorption interaction. 

Because neutrino interactions are expensive to compute on-the-fly we 
have implemented interpolation tables. Neutrino interactions depend on 
the neutrino momentum components and the state of the fluid with which 
the neutrinos interact. A grid of neutrino energies and angles is fixed, but 
the fluid density n, temperature T, and electron fraction Y e vary throughout 
the simulation. Hence we employ tables that may be interpolated in n, T, 
and Y e . Particularly for neutrino scattering and pair interactions — which 
depend on neutrino states before and after the collision — the interaction 
kernels arc of high dimensionality, requiring a globally distributed table. 
On each processor a local table is constructed, which contains a copy of 
each n, T, and Y e vertex required by the zones for which that processor is 
responsible. As n, T, and Y e in a processor's zones evolve, the relevant ver- 
tices are pulled from the global table as needed. For each zone we construct 
a "cube" of pointers to the eight vertices surrounding the zone's values of 
n, T, and Y e . This cube is then used for the necessary interpolations. 
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7. Outlook 

We have made a promising start on GenASiS, a new code being developed 
to study the explosion mechanism of core-collapse supernovae. Our plan 
is to include all the relevant physics — including magnetohydrodynamics, 
gravity and energy- and angle-dependent neutrino transport — in a code 
with adaptive mesh refinement in two and three spatial dimensions. Paral- 
lelization and implementation on the adaptive mesh are not yet complete, 
and the physics components have not yet been fully integrated; but steady 
progress and the successful completion of test problems give us confidence 
that we are well on our way towards a tool that will provide important 
insights into the supernova explosion mechanism. 
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